Nonlinear dynamo action in a precessing cylindrical container 
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It is numerically demonstrated by means of a magnetohydrodynamics (MHD) code that precession 
can trigger the dynamo effect in a cylindrical container. This result adds credit to the hypothesis that 
precession can be strong enough to be one of the sources of the dynamo action in some astrophysical 
bodies. 
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Simulating numerically or reproducing experimentally 
the dynamo action in astrophysical bodies is challeng- 
ing in many respects and the question of converting ki- 
netic energy into magnetic energy through a realistic 
forcing mechanism is seldom addressed either numeri- 
cally or experimentally. For instance, numerical simu- 
lations are often restricted to kinematic dynamos where 
fluid flows are prescribed for simplicity reasons and bear 
little resemblance to reality. The three fluid dynamo 
experiments [U |8l |9j that have been successful so far 
do not model astrophysical dynamos either, since the 
external pumps or internal impellers which have been 
used in these experiments have no astrophysical coun- 
terparts. As far as natural forcing mechanisms are con- 
cerned, buoyancy and precession are believed to be pos- 
sible sources of energy for the geomagnetic dynamo. The 
precession hypothesis has been formulated for the first 
time in [7J and has since then been actively investi- 
gated from the theoretical, experimental and numeri- 
cal perspectives t2j |5l l6j. To the best of our knowl- 
edge, however, it seems that it is only recently that nu- 
merical examples of precession dynamos have been re- 
ported in spheres [W, llj and in spheroidal cavities [12j. 
While spheres and spheroids are relevant for planetary 
dynamos, cylindrical containers seem more convenient for 
experimental purpose. In this respect we have in mind 
the large scale MHD facility DRESDYN currently being 
built at Helmholtz-Zentrum Dresden- Rossendorf in Ger- 
many where, among other things, the action of precession 
will be tested on cylinders (F. Stefani, personal commu- 
nication). The objective of the present Letter is to report 
numerical evidences supporting the idea that precession 
is indeed a potent mechanism to drive dynamo action in 
cylindrical containers. 

The conducting domain considered in this letter is a 
cylindrical vessel C of radius R and length L. The vessel 
contains a conducting fluid and is embedded in vacuum. 
The solid walls of the vessel are assumed to be so thin 
that their influence is henceforth neglected. The con- 
tainer rotates about its axis of symmetry with angular 



velocity Qr^^ and is assumed to precess about a second 
axis spanned by the unit vector forming an angle a 
with e^, (0 < a < tt). The precession angular velocity is 
QpGp. A cylindrical coordinate system about the axis of 
the cylinder is defined as follows: the origin of the coor- 
dinate system is the center of mass of the cylinder, say 
O; the Oz axis is the line passing through O and parallel 
to e^; the origin of the angular coordinate ^ (0 < ^ < tt) 
is the half plane passing through O, spanned by and 
Gp, and containing rtpSp. The third coordinate, denoted 
r, is the distance to the Oz axis. 

We denote hy C = R and U = RQr the reference 
length and velocity scales, respectively. The fluid den- 
sity, p, is assumed to be constant and the reference pres- 
sure scale is V := pU'^ . The magnetic permeability is 
uniform throughout the entire space, /io, and the elec- 
tric conductivity of the conducting fluid is constant, ctq. 
The quantities /io and <Jo are used as reference magnetic 
permeability and electric conductivity, respectively. The 
reference scale for the magnetic field is chosen so that 
the reference Alfven speed is 1, i.e., 
are left with five non-dimensional parameters: one geo- 
metrical parameter LjR (aspect ratio); two forcing pa- 
rameters OL (precession angle) and e = Q^pjO^r (precession 
rate); and two fluid parameters, namely the Ekman num- 
ber E = u/R'^Qr (where u is the kinematic viscosity) and 
the magnetic Prandtl number Pm = vijlqctq. We finally 
define the kinetic Reynolds number Re = 1/E and the 
magnetic Reynolds number Rm = PmRe. 

The non-dimensional set of equations that we consider 
is written as follows in the precessing frame of reference: 

5^u+ (u-V)u + 2£:epXu + Vp = ^Au + f, 

Re 

Vu = 0, 

dth-Vx(u X h) = ^Ah, 
Rm 

Vh = 0, 

where u, and h are the velocity field, the pressure, 
and the magnetic field, respectively. In the following 
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we consider three different modes to solve these equa- 
tions: (i) The incompressible Navier-Stokes mode; (ii) 
The Maxwell or kinematic dynamo mode; (iii) The non- 
linear magnetohydrodynamics mode (MHD). In Navier- 
Stokes mode the source term f is set to zero and h is 
not computed. In Maxwell mode, only the induction 
equation is solved assuming that the velocity field u is 
given. In MHD mode the full set of equations is solved 
and the source term f is the Lorentz force per unit mass, 
f := (Vxh)xh. The no-slip boundary condition on the 
velocity field is written as follows in the precessing frame 
of reference: u = e^i at r = 1 and u = tgq at z = ±1. The 
magnetic field is represented as the gradient of a scalar 
potential in the vacuum, V^. The magnetic boundary 
transmission conditions enforce that the magnetic field 
is continuous across the walls of the vessel, say S, i.e., 
h|s = V(/)|s. 

The above equations are solved numerically by means 
of a code which is specialized to axisymmetric domains 
and has been presented in details in [3, 4J. The code 
is called SFEMaNS for Spectral/Finite Elements for 
Maxwell and Navier-Stokes equations. It is an hybrid 
algorithm that uses finite element representations in the 
meridian section of the axisymmetric domain and Fourier 
representations in the azimuthal direction. The magnetic 
field is represented as a vector field in the conducting re- 
gion and as the gradient of a scalar potential in the insu- 
lating region. SFEMaNS can account for discontinuous 
distributions of electric conductivity or magnetic perme- 
ability and all the required continuity conditions across 
the interfaces are enforced using an interior penalty tech- 
nique. The solution technique is parallel and paralleliza- 
tion is done with respect to the Fourier modes. 

The typical spatial resolution in the meridional plane 
of the conducting domain is Ax = 1/160. We take 32 
Fourier modes (m = 0, ...,31) for Navier-Stokes runs 
and 64 Fourier modes (m = 0, . . . , 63) for MHD runs. 
The typical time-step is At = 0.001. The grid is non- 
uniform in the vacuum with Ax = 1/160 at the cylinder 
walls and Ax = 1 at the outer boundary of the numeri- 
cal domain, which is a sphere of radius ten times larger 
than that of the cylinder. A typical MHD run requires 
about 1000 CPU hours per rotation on 64 processors on 
an IBM-SP6. 

Let us first briefiy recall what is observed in a typical 
precessing fiuid experiment starting with the fiuid at rest, 
(see e.g. [5^, 'S]). The vessel is first set in rotation with- 
out precession. The fiuid motion is then governed by the 
formation of a viscous Ekman boundary layer during the 
acceleration ramp. The resulting fiow is a stable solid 
rotation independently of the strength of the accelera- 
tion phase. Once precession is applied, the Coriolis force 
generates an axial motion of the fiow supported by the 
Fourier mode m = 1. When Re is large enough, the fiow 
undergoes a transition from laminar to turbulent even 
for small precession rates and small angles [5j . The range 



£ G [0.1, 0.15] has been shown in [H] to maximize the ratio 
of axial to transverse energy in a cylinder of aspect ratio 
2 in the range Rm G [500, 5000] when a = 7r/2. Although 
a detailed study of the various transitions between these 
hydrodynamic regimes is interesting per se, due to lim- 
ited numerical resource we reduce the dimensionality of 
the parametric space to one aspect ratio, L/R = 2, one 
precession angle, a = 7r/2, one precession rate, e = 0.15, 
and only two values of Re G {1000, 1200} and four values 
of Rm e {600, 800, 1200, 2400}. 

We start our investigations with a Navier-Stokes run 
at Re = 1000. The initial velocity field is the solid ro- 
tation in the precessing frame: uq = e^xr. The onset 
of the axial circulation induced by precession is moni- 
tored by recording the time evolution of the normalized 
total kinetic energy K{t) = \ /^u2(r,t)dr/Ko and nor- 
malized axial kinetic energy Kz{t) = \ f^u1{r,t) dr/Ko 
where Kq = ^ f^'^o dr is the kinetic energy of the ini- 
tial motion. The time evolution of K{t) and Kz(t) for 
t e [0, 272] is reported in fig. [l] The time t = 272 cor- 
responds to 43.3 rotation periods. After a transient that 
lasts 5 rotation periods and peaks at two rotation pe- 
riods, the axial kinetic energy reaches a plateau value 
Kz ~ 0.1. Meanwhile, the total kinetic energy decreases 
and reaches a plateau value K ^ 0.42 after 5 rotation 
periods also. To enrich the dynamics of the system we 
have restarted the computation at t = 72 (i.e., 11.5 ro- 
tation periods) and increased the Reynolds number to 
Re = 1200. The time evolution of K{t) and Kz{t) for 
t e [72,342] and Re = 1200 is also reported in fig. [l] 
The time evolution of the total kinetic energy shown in 
fig. |l(c)| presents doubly periodic oscillations with one 
long period of about 8 rotation periods and one small 
period of about one rotation period. The short period 
oscillations correspond to energy exchanges between the 
north and south halves of the container, with a period of 
2 rotation periods. The energy exchange mechanism is 
visible in fig. 1(d) [ where we have reported the time evo- 
lution of the kinetic energy of the north and south halves 
of the cylinder for t G [312,342]. Similar oscillations be- 
tween north and south hemispheres have been reported 
to occur in a spheroidal cavity in [12j. 

We now solve the full MHD system using as initial ve- 
locity field the velocity computed at t = 192 during the 
Navier-Stokes run at Re = 1200. The initial magnetic 
field is defined as follows to trigger efficiently the dy- 
namo instability. For the Fourier modes m G {0, 1}, the 
homogeneous Dirichlet boundary condition on the scalar 
potential is replaced hy (j) = 0.05zf{t) for m = and 

(j) = 0.05r/(t) for m = 1, where f{t) = (l 



with 



0^ and 



2 • 



For m > 2, the amplitude 



of each Fourier mode of the initial magnetic field com- 
ponents is set to 10~^. Various MHD runs are done at 
Re = 1200 for different values of the magnetic Prandtl 
number. The onset of dynamo action is monitored by 
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(a) total kinetic energy K, 
Re = 1000, Re = 1200 
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(b) axial kinetic energy Kz, 
Re = 1000, Re = 1200 
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(c) zoom of (a), f G [50,342], (d) zoom of K/2 in (a) and north 
Re = 1000, Re = 1200 and south kinetic energies. 

Re = 1200, t e [230,342] 

Figure 1: Time evolution of the total kinetic energy K, axial 
kinetic energy Kz, and total north and south kinetic energies 
as indicated. 



recording the tirae evolution of the magnetic energy in 
the conducting fluid, M{t) = ^ f^h'^{r,t)dr/Ko. Dy- 
namo action occurs when M{t) is an increasing function 
of time for large times. The time evolution of M for 
Pm G {2, 1, |, ^} are shown in fig. 2(a) 



The runs at 

Pm G {1, |, ^} are done by using the velocity and mag- 
netic fields obtained from the run Pm = 2att = 211as 
initial velocity and magnetic fields. The flow is above dy- 
namo threshold for Pm = 1 and Pm 
cal for Pm - 

gives the critical magnetic Prandtl number Pm* ^ 0.625 
corresponding to the critical magnetic Reynolds number 
i?m* ^ 750. 
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Figure 2: Time evolution of the magnetic energy M in the 
conducting fluid (a) in the linear regime from t — 192 for 
Re — 1200 and various Rm as indicated (in lin-log scale) 
and (b) in the nonlinear regime from t — 192 to t = 241 
{Re 1200, Rm 2400), from t = 241 to t = 271 
{Re 1200, Rm 1200) and from t = 271 to t = 292 
{Re^ 1200, Rm^mO). 



We now wish to observe the nonlinear saturation and 
evaluate the impact of the magnetic Prandtl number on 
the nonlinear regime. To reach nonlinear saturation in 
reasonable CPU time, we have used as initial data for 
the velocity and magnetic fields the velocity and mag- 
netic fields from the IvlHD run Pm = 2 at t = 217. 
The velocity field has been kept unchanged but we have 
multiplied by 300 the amplitude of the Fourier modes 
m = 0, . . . , 5 of the magnetic field. The time evolution 
of the magnetic energy of this nonlinear run in the time 
interval t G [192, 241] is shown in fig. |2(b) We observe 



that M grows smoothly until t = 222 and begins to os- 
cillate thereafter. The ratio M/K is observed to be of 
order 10~^ during the nonlinear oscillating regime. After 
restarting the IVIHD run at t = 241 with Pm — 1 and 
running it until t = 271, we observe that the dynamo is 
still active. After restarting the IvIHD run at t = 271 
with Pm — \ and running it until t = 292, we observe 
that the dynamo dies in a short time lapse, indicating 
that the dynamo bifurcation is not sub-critical for this 
set of control parameters. This experiment confirms the 
interval \ < Pm^ < | for the critical magnetic Prandtl 
number for dynamo action which has already been ob- 
served in the linear regime. 

Tilgner [lOj has observed that unsteadiness and break- 
ing of the centro-symmetry of the flow facilitate dynamo 
action. A similar observation has been made in [12], and 
dynamo action is reported therein to occur when cyclic 
oscillations of the kinetic energy between the north and 
south halves of the spheroidal cavity occur. Although 
the loss of centro-symmetry is not a necessary condition 
for dynamo action, we now want to test this idea in the 
present cylindrical setting. 

The loss of centro-symmetry of the velocity field can 
be monitored by inspecting its symmetric and antisym- 
metric components: Us{r,t) = ^(u(r, t) — u(— r, t)) and 
Utt(r, t) = ^(u(r, t) + u(— r, t)). In the Navier- Stokes sim- 
ulations reported below, we monitor the loss of centro- 
symmetry by inspecting the time evolution of the asym- 
metric kinetic energy Ka{t) = \ u^(r, t)(ir/i^o and we 
define the asymmetry ratio ra{t) = Ka{t) / K{t). The 
computations reported below have been done on centro- 
symmetric grids, but centro-symmetry is not otherwise 
enforced. 

The time evolution of the asymmetry ratio is shown 
in fig. |3]for the precessing cylinder at Re = 1200 in the 
time range t e [72, 342] (dotted line). The ratio varies 
between 0.004 and 0.01 when the nonlinear regime is well 
established, i.e., t > 220. In order to evaluate the impact 
of the dynamo on the symmetry of the flow, we have 
started the MHD run at t = 192 with Pm = 2 (i.e., 
Rm = 2400). The time evolution of is shown in solid 
line in fig. [3] Note that the solid and dotted lines coin- 
cide since the dynamo regime is linear in the time interval 
t G [192, 217] and the magnetic field is too weak to have 
an impact on the energy ratio r^. We have restarted the 
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MHD run at t = 217 after raultiplying the amplitude of 
the magnetic field by 300 as already mentioned. The ra- 
tio ra clearly departs from its Navier-Stokes value there- 
after as seen in the figure. At saturation, oscillates 
between 0.09 and 0.011; these values are slightly greater 
than those reported in [101 ^ precessing sphere. We 
have restarted the MHD run again at t = 241 after reduc- 
ing the value of Pm to 1, thereby reducing the magnetic 
Reynolds number to Rm = 1200. The asymmetry factor 
is not dramatically affected by the change, as seen on the 
figure. We have finally restarted the MHD run at t = 271 
after reducing the value of the magnetic Prandtl number 
to ^. As expected the dynamo dies and decreases to 
its hydrodynamical level. 
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Figure 3: Time evolution of the asymmetry ratio ra at 
Re = 1200 for t G [72, 342] in Navier-Stokes regime and Re = 
1200, Rm = 2400 for t = [192,241], Re = 1200, Rm = 1200 
for t = [241, 271], and Re = 1200, Rm = 600 for t = [271, 292] 
in MHD regime. 
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Figure 4: Time evolution of the magnetic energy M at 
Re = 1200 and Rm = 1200 for t G [211,229] in MHD mode 
(denoted as 'MHD'), in Maxwell mode with the velocity frozen 
at t = 211 (denoted as 'MAXWELL') and in Maxwell mode 
with the symmetrized velocity frozen at t = 211 (denoted as 
'MAXWELL SYM'). 



m = 1, 2, 3. 




(a) from the side (b) from the top 



In order to study the impact of the centro-symmetry 
and the unsteadiness of the flow on the dynamo action, 
we have performed two Maxwell runs at Rm — 1200 
with the following characteristics: (i) the velocity field 
at Re — 1200 is frozen at t = 211, (ii) the velocity field 
at Re = 1200 is frozen at t = 211 but only its symmetric 
component is retained so that the resulting velocity field 
is centro-symmetric. The time evolution of the magnetic 
energy of the MHD run and the two Maxwell runs (i) 
and (ii) are shown in fig. [4j It is remarkable that, in the 
two considered kinematic runs, the dynamo keeps grow- 
ing with a rate similar to that of the MHD run. These 
computations show that neither the temporal oscillations 
nor the flow asymmetry play a crucial role on the dynamo 
action in the precessing cylinder at Rm = 1200. 

A snapshot of the vorticity and magnetic lines at 
Re = 1200, Rm = 2400 is shown in fig. (Sj We ob- 
serve a central S-shaped vortex which is deformed by 
the precession and reconnects at the walls through vis- 



cous boundary layers, (see fig. 5(a)). The magnetic field 
lines exhibit a quadrupolar shape which is best seen in 



the vacuum from the top of the cylinder (see fig. 5(b)). 



The magnetic energy is dominated by azimuthal modes 



Figure 5: Snapshot at t = 241 for Re = 1200, Rm = 2400 
showing vorticity field lines (red) and magnetic field lines col- 
ored by the axial component in the cylinder (grey/black for 
positive/negative hz component). 

Forty years after the promising experiments with liquid 
sodium by Gans [2j, we have numerically demonstrated 
dynamo action in a precessing cylindrical tank. The bi- 
furcations through symmetry breaking and cyclic time 
dependence are similar to those already observed in dy- 
namo flows in spherical or spheroidal precession driven 
cavities. There is however a large gap between the con- 
trol parameters used in the present simulations and those 
achieved in experimental set-ups and planetary dynamos, 
where E = 1/Re and Pm are many orders of magni- 
tude smaller. Following this preliminary evidence for 
dynamo action, two further steps appear now as most 
urgent: (1) studying parity breaking and unsteadiness 
by varying the forcing parameters (precession angle and 
rate); (2) searching for a scaling law for the critical mag- 
netic Reynolds number as a function of the hydrodynamic 
Reynolds number. Such a relation has been proposed by 
Tilgner in a precessing sphere [lOj, who argues that it 
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is the asymmetric part of the flow that plays a key role 
in the dynamo. The research program (2) will be time 
consuming as it will necessitate large scale computations 
to explore a wide range of Reynolds numbers. It will also 
require to develop nonlinear stabilization techniques to 
simulate small scale viscous dissipation. A major step in 
the understanding of precession dynamo will hopefully 
be achieved in the near future with the construction of 
the large scale MHD facility DRESDYN at Helmholtz- 
Zentrum Dresden-Rossendorf (Germany). The coopera- 
tion between simulations and experiments will lead to a 
better understanding of natural dynamos, including the 
geo dynamo. 
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